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Abstract 

We use a systematic method which allows us to identify a class 
of exact solutions of the Flierl-Petvishvili equation. The solutions are 
periodic and have one dimensional geometry. We examine the physical 
properties and find that these structures can have a significant effect 
on the zonal flow generation. 
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1 Introduction 

The problem of coherent structures has been extensively investigated in con- 
nection with the drift wave modes in tokamak. Numerical simulations and 
experimental observations have provided strong evidence of the presence of 
long-lived, cuasi-coherent structures even in deep turbulent regimes. The 
theoretical models (close to similar descriptions in the physics of fluids, at- 
mosphere and ocean) emphasize the role of two types of nonlinearity : scalar 
(or Korteweg deVries-type) and vectorial (convection of the vorticity) . Both 
are able to support vortical flow in a form of coherent and long lived struc- 
tures. The problems of generation from initial conditions and the stability 
of monopolar or multipolar vortices are subjects of intense research and one 
of the notable result is that not all of these structures can be expected to 
be solitons in the sense of the Inverse Scattering Method. Usually they are 
called solitary waves or solitary vortices. A review, mainly oriented to plasma 
physics applications, has been done by Horton and Hasegawa pQ. 
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The ion drift wave equation has a distinct dynamical character according 
to the space-time scales involved. For shorter scales the Hasegawa-Mima- 
Charney equation is obtained and the dynamics exhibits dipolar structures 
on the Larmor radius scale. On larger scales the scalar (or KdV-type) nonlin- 
earity is prevailing and the structures are monopolar. Both are not solitonic 
but very robust and long lived. In the latter case the one- dimensional version 
of the equation can be reduced to the modified Korteweg De Vries (mKdV) 
equation |2j and has been derived in various contexts in plasma physics. 
Tasso [Sj and later Petviashvili [4 j have derived versions of this equations ap- 
plicable to the case where there is a strong temperature gradient. It has been 
shown [H] that the equation proposed was in fact exclusively dependent 
of the density gradient. Later the equation has been rederived along with 
a careful analysis of the scales involved [7], resolving a controversy on the 
role of the temperature gradient. In the study of ocean flows Flierl |Hj has 
independently formulated an equation with the same structure. A one dimen- 
sional version of the equation has been solved by Lakhin et al. on an infinite 
domain [S], obtaining as solution the KdV soliton. The two dimensional case 
has been examined by Kadomtsev and Petviashvili [5| and by Petviashvili 
et al. POj using trial functions to extremize a functional derived from the 
equation. The solution they found is vortical monopolar. Boyd and Tan [llj 
have found a monopolar solution expressed as a sum of 45 terms depending 
on the radial coordinate r through trigonometric functions. They have also 
proved the non-existence of stable non-axisymmetric solutions to the station- 
ary Flierl-Petviashvili equation ^2j. There is a strong connection between 
the stationary Flierl-Petviashvili equation and the Zakharov-Kuznetsov equa- 
tion PHI, a two-dimensional generalization of the KdV equation. This has 
also been investigated in Refs. [T^j and jT^j where multidimensional nonlin- 
ear wave structures (of the electrostatic drift wave branch) in inhomogeneous 
plasmas were studied by combining the eigenvalue problem in inhomogeneity 
and the 2D nonlinear vortex equation. For a two-dimensional generalization 
of the KdV equation numerical results are available For shorter refer- 
ence it is used also the name Regularized Long Wave equation, RLW, since 
the equation proposed by Peregrine J7j can be brought to the same form. 

The stationary two-dimensional equation derived in the studies mentioned 
above has a very simple structure, and a one-dimensional solution is known, 
the KdV soliton. We will examine the possible extension to an analytical 
closed form of a two dimensional solution. We are motivated by the studies of 
plasma ion mode instabilities, in particular the generation of radially localised 
layers of sheared flow (zonal flow). 
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2 Derivation of the equation 



The ion instabilities are dominated by the nonlinearity related to the polar- 
ization drift. When the temperature gradient is very small the quasi-three di- 
mensional geometry assumed in the derivation of the Haswgawa-Mima equa- 
tion leads to suppression of the classical convection nonlinearity although the 
presence of the later can still be important in a fluctuation model. The ion 
polarisation drift nonlinearity may support vortical flows of which certain 
states are stable and coherent (actually not solitonic). The derivation of a 
nonlinear equation in this case starts by the ion density continuity 

d \ 

^ + v i± - V ± Jn i + n (V ± -v) = (1) 

and assuming neutrality and adiabatic response of the particles 

rii « n e « n — 

J- e 

The ion velocity in two dimension is 

' --V ± )V ± cp (2) 



B VtiB \dt B 

The Eqs.(JH) and (J2J) lead to a nonlinear equation for the electrostatic po- 
tential. However, in these two equations there is still much freedom and a 
detailed analysis of the time and space scales leads to different particular 
forms of this equation. This is a multiple space and time scale analysis 

U = e l t 



for £ < 1. It is possible to change to a moving frame, by introducing a 
translation velocity, u. 

r] = y — ut 

The speed u can also be a function on certain space and time scales. 

We recall the main steps of the analysis performed in Ref. j7j. In the 
scaling 

if = Etp 1 + e 2 <p 3 + ... 

with 

ip x = (fx (x , X U X 2 , T] , r} x ,7} 2 , t , ti, t 2 , ...) 
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u ~ ui = O (e) 



n = n {x 1 ,x 2 , ...) 
T e = T e (x u x 2 ,...) 

In this scaling no and T e have variations on only "large" space scales, and 
the potential (which is of e -small amplitude) varies on the small space scale, 
xq and slow time scale, t an d, of course, on higher scales. The characteristic 
space scale for the potential is then p s which is the spatial extension of the 
dipolar vortex. The normalisation of the parameters can be done according 
to these scales 



, &p T e v i± 

<p = — , 1 — — , v — — 
T T c s 

x' = x/p s , 7]' = t]/ Ps , t' = n~H 



and dropping the primes, one obtains 

1 _ 



_d_ 



T( Xl ) 
d 



+uij—V ±0 (l) 1 



Ui 



T( Xl ) 



+ 



[(-V ±o 0i x n) • V 



drj 



where 

= ! _ 1 dn 
Kn ~ K ~^ ^x~ 

is varying on the larger space scale, x±. 

This is the Hasegawa-Mima equation governing the space variation of the 
potential <p 1 on the smaller space scale, xq and i] of the order p s and on times 
of the scale t±. On these scales several simplifications are obvious, since the 
temperature and the density gradient lengths can be taken constants 

K n / U i ~ constant on x 
T ~ constant on xq 

A different dynamical equation is obtained on other space-time scale 

<jy = e 2 <f) 2 + e 3 (j) 3 + ... 
02 = 02 (xi,x 2 , ...,r] 1 ,V2, ■■■,h,t 6 , ...) 
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n Q = n (x 2 ,x 3 , ...) 
u = u 2 ~ O (e 2 ) 
T = T(x 2 ,x 3 ,...) 
The following combination has to be of a certain order in e 

T (*„) u 2 U {£ > 
Then the equation on the time scale of order 5 is 

9 02 9 _ 2 , 



9t 5 T(a; 2 ) 



-«2 



1 . «n (^2) 



T (x 2 ) U 2 

d(f) 2 



d(f> 2 
drix 



+ K T {X 2 )(p 2 -^— 

= 

The operator of Laplacean in two dimension acts on the larger space scales 

^2 d 2 d 2 
VI = 1 

dx\ dr/f 

The dominant space variation is here on the scale 

Ps 

xi = — 

e 

which is much larger than the dipolar vortex scale. The condition imposed 
by this ordering is 

dn n 1 

k t = « 

OX\ u 2 

In this range the nonlinear equation is dominated by the scalar nonlin- 
earity. This equation leads, after one integration over the r/ 1 coordinate, to 
the Flierl-Petviashvili equation. 

Spatschek et al. [Zj also discuss intermediate scalings, where both the 
scalar and the vectorial nonlinearities are present. The stability analysis 
implies the concept of structural stability, where the dynamics governed by 
one of the equations obtained above (Hasegawa-Mima or scalar nonlinear) 
is perturbed with a term that is of the other type. It has been shown that 
the dipolar vortices are broken into separate monopolar vortices, while the 
monopolar vortices are structurally stable. To these consideration one should 
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also add the numerical results from the collision of monopolar vortices of the 
Flierl-Petviashvili equation (or Zakharov-Kusnetsov) [TT], [!§]• The 

monopolar vortices are stable and are not destroyed by collisions, however 
the form and the amplitides are perturbed, showing again that they are not 
exact solitons. 



3 The scalar nonlinearity equation 

In plasma physics applications, the Flierl-Petviashvili equation has the form 

A0 = acj) - (3(j) 2 (3) 

where a and (3 are physical parameters, i.e. functions depending on (x,y). 
This equation is obtained as the stationary version of the equation which has 
been expressed in the system of reference moving with the velocity u 

y^y-ut 

In the following, the coordinates x and y are the coordinates in the moving 
system. It is assumed that y is the poloidal direction and x is the radial 
direction in tokamak. 

In the tokamak context the equation has been frequently simplified 

d 2 d 2 



dx 2 dy 

leading to 



° 2 ° a<P-(3<p 2 (4) 



dy 2 

whose solution (Lakhin et a/.[S]) is 

cosh ( ^y ) 

The dependence on x is only parametric here, via the coefficients a and j3. 
This has the same form as the KdV soliton solution. 

Assuming that a and (3 are constants, this solution can be extended 
indefinitely in the transversal (x) direction, as a ridge propagating in the y 
direction whose section is given by Eq.(J3J. If the plasma is homogeneous, 
this ridge can be rotated arbitrarly in plane. We have then a family of quasi 
two dimensional solutions to the Eq.(j3J). 

We will discuss the possibility to find other solutions of the two-dimensional 
version of the equation. 
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4 Starting from the one dimensional problem 



We can directly integrate the one dimensional version of the equation (j3J). 
Actually, the solution (jHJ) is one of the possible solutions arising from the 
direct integration. For a reason that will become clear later, we write the 
one dimensional version of the equation (jSJ) in the form 

d 2 6 



d (irj) 2 



am — 



where it] can be y as in Eq.(j3J). Multiplying by d<p/dr) and integrating once 
we have 

1 f d(f)\ 2 a 2 fl 3 
-2UJ = 2 +K 
where k is a constant. Then we have the integral 



Jf <f - a<p 2 - 2k 

The integral is elliptic and its inverse <fi (77) has a closed analytical expression 
for all k. Particular expressions can be written according to the reality of 
the roots of the third degree polynomial under the square root. For k such 
that all roots (a > b > c) are real, we have, for (ft > a 

/ —j^^^^^=^^^^^^ — 9 sn (sin 93, k) 

U V - a ) -b)((f)-c) 

= gF( v ,k) 

where sn is the Jacobi elliptic sinus, F is the incomplete elliptic integral of 
the first kind and the notations are 



if = arcsm 
2 



-y/a — c 
b-c 



9 

k 2 

a — c 

This expression can be easily inverted to obtain as function of 77. An 
analoguous result is obtained when k is such that the root a is real and b and 
c are complex. The integral is written, for > a 

= = gen -1 (cos (p, k) 

= gF(<p,k) 
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where 



(6 + 6) (6-6) 2 
61 " —2~> ^ ~ r~ 

g = — , A 2 = (61 - a) 2 + a 2 
v A 



2 A + 61 - a 

fc " 2A 

cf) — a — A 

09 = arccos 

(f) - a + A 

Deatils on these Jacobian elliptic functions can be found in Ref.|18j. 



5 Constructing the solution from singulari- 
ties 

In view of the application to the description of plasma drift wave, a two- 
dimensional solution would be very useful. We already dispose of the one- 
parameter family of solutions with geometry of fronts, which has been ob- 
tained by simply translating a one-dimensional solution along the transversal 
direction. We examine a class of solutions that can be seen as an extension 
from the one-dimensional ones. We propose a more systematic approach of 
construction that makes more explicit the nature of the restriction leading 
to this quasi-one dimensional geometry. 

We will start from the one dimensional model, taking the coefficients 
constants. The only analytical formula for a solution of (JHJ) is |Hj 

(y) = 0o sech2 (iv) 

where 

3a 



7 = 2 

(note that in Lakhin et al. the coefficient must be corrected by dividing with 
2). We will use the following relation 



sech 2 (z) = — cosech 2 I z 



ITT 
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and the expression of the cosech function as a series using its pole singulari- 
ties. Then we obtain 

X j 

(y) = o sech 2 (jy) = -</> V (6) 

,=_«, {jy - f + zZvr) 

We are looking for solutions (x, of the 2D equation and we try to 
build them from the motion of the pole singularities of the one- dimensional 
solution. More specifically, we will make an ansatz for the form of the 2D 
solution based on a particular choice of poles in the complex y plane, as 
suggested by Eq.©. We will assume that the poles have positions that 
depend on the other coordinate, x. Imposing that the function constructed 
in this way verifies the equation we obtain a set of differential equations and 
constraint conditions for the positions of the poles in the complex y plane. 
The poles evolve with x as time-like variable. Expressed in other terms, we 
look for a function that is meromorphic in the complex y plane and whose 
poles depend on x. We find, following other similar approaches, that the 
solution is an elliptic function, i.e. a doubly periodic meromorphic function 
of y, for all values of x. 

This procedure has been deveoped in the context of the exactly integrable 
differential equations like Kortweg de Vries. We follow closely the methods 
exposed in papers of Choodnovsky [20] , Thickstun [21] . The work of Decon- 
inck and Segur [22] is particularly relevant for our problem. 

The following ansatz is suggested by the theory of r-functions in the 
integrable equations context jTHj 

d 2 

4>{x,y) = 2— \nT[x,y) (7) 



where 



(x, y) = c TT ( 1 - — ) exp ( — ] 

fe = i V VkJ \VkJ 



On the other hand, there is also the suggestion from Eq.© that the function 
<p is periodic on the imaginary y axis. We choose a periodicity iD (compare 
with in in the above equation) and take iV poles in each domain of periodicity. 
Then r is 

N oo j 

We insert this expression in Eq.fJJJ and take care of the constants 

N oo 

(x, y) = -20 o E T~( , - m1 2 ( 10 ) 



=i i- 



~1 [7 (y - Vn) - UD]' 
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6 The dynamics of the singularities 



6.1 Extending by inclusion of the x coordinate 

The x dependence in this equation comes from the dependence on the variable 
x of the position y n (x) of the poles in the complex y plane. Now we impose 
that this form of 4> verifies the equation. 



dx 2 dy 2 

N 



JV OO y 

-47*0 E E ( 

n=l /=— oo ^ 



dy n (x) 



dx S±t V dx ) h(y-y n )-UD} 



d2( t> A ^ ^ f 1 



2 



1 



n=l l=-oo 

3 ( dyn (x) 

V dx J [ 7 (y - y B ) - iW] 1 
o , N oo 

70o ^^[7(y-yn)-^] 3 

(9 2 iV oo ^ 

We replace these formulas and those for and 2 in the equation. We 
will examine the neighborhood of one of the poles, taking 

y = y P (x)+e 

where e is a small quantity and p is one of the N poles . Expanding all terms 
in the equation in e we equal to zero the coefficients of the same powers of 
e. This will give us the equations we have to impose to y n (x). 
We have, with e — > 0, 

, (~20o) 1 

= —^T 2 + 

N oo 



+ (- 2 <^o) E + 1- 2 ^ E E 



Z^O n^P 
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2 



3V 

dx 2 



v-40o7) 



d 2 y p (x) \ 1 



+ (-4</»o7) 



dx 2 J 7 3 5 3 
d 2 y p (a 



+ 



OO 



1^0 



+ 



n^p 

+ (-i2^)f (%M) 2 f 

n=l ^ 



'd 2 y n (x) 



+ 



n^p 



1= — oo 



[7 (y P - Vn) - ud] 



d 2 



dy 2 



(" 12 ^o7 2 ) ^ + 



i=— 00 

N 



+ (-12^) £ £ 



n=l i=-oo 



[7 (Vp ~ Vn) ~ UD] 



2 = 



+ 



7% 2 



N 00 



2 y -^+2V y 



n=1 [7 {Vp-Vn) -UDY 

n^p 



+ 



+ 



N 00 



1 



=—00 
L 1^0 



- ilD ) t^l [7 (y P - Vn) - iwy 

njtp 



We have to collect the terms containing the same powers of e, begining with 
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the highest. For easy identification, we show separately the contributions of 
each term of the original equation. 
For 1/e 4 : 



(-120o7 2 ) 



dy p (x) 



dx 



7 % 4 



+ 



1 



(-120o7 2 ) ^ 4 
(-a) x + 





This gives the equation 



dy P (x) \ 2 = 1 o /3 _ i 
dx 



or 



% (a) 



3 7 2 

2 



For 



(-40o7) 
+ 

(-a) x + 
/3 x 




d 2 y p (x) \ 1 



dx 2 J 7% 



+ 



The resulting equation is 



d 2 y P (x) 
dx 2 



Finally, the coefficient of 1/e 2 : 
+ 
+ 



(-a) x 



-20 o ) 1 

7 2 e 2 



+ 



1 







AT oo 



(11) 

(12) 



(13) 



2 v -^+2V y f 



2^0 
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This results in the following constraint 

oo j N oo 

1 

~3 

From the Eqs. (fT3|) and (|T2"J) it results that the "trajectories" y p (x) are 
linear on x. We have 

%M = ±1 (is) 
ax 

or 

y p (x) = ±x + c p (16) 

where constants. 



6.2 Discussion of the constraint equation 

We have obtained the equations and the constraint using an expansion around 
a singularity y p . The choice is arbitrary and we can repeat the calculations 
using another singularity. Obviously, the form of the equations will not be 
changed. Finally, we will obtain a number of constraints equal to the number 
of singularities. 

A simpler form of the constraint can be written. 



V - = - 2 V - 

^{-iWf (-iD) 2 

1^0 



We introduce the notation 



and we have 



= !2 ( _ ^ 

z pn — jj \Vp Vnj 



1 7T 1 



n2 



n^p 



The second sum can be written 



1 1 



13 



and the constraints become, for p = 1, N 



N 

E 

n=l 



7 2 (Vp - J/n)' 



/ / Up Un 

iDpf 



Vn ~ Up 



ti 2 -D 2 



In these formulas, ip' is the first derivative of the Euler ^si-function. This 
form can be useful if one wants to examine numerically the validity of a 
particular choice of constants c p . 

Alternatively we can use the expansion for the square of the cosech func- 
tion, as before. 



E 



7T 



b{yp-y n )-UDf D 2 



= -—cosech 



l=— oo 

and the form of the constraint becomes 

^2 N 

^2 cosech 2 (y p - y n ) 



7T7 

~D 



(Vp - Vn) 



D 2 



n=l 
n^p 



The following identity exists 

s-1 



71 



3D 2 3 



(17) 



E 

k=i 



9 / kn 
cosec — 

V s 



s 2 -l 



(18) 



Comparing our equation with the identity and using the elementary relation 

cosech 2 (ix) = —cosec 2 (x) 
it is suggested to identify 



D 



71 

l(Vp-Vn)/i = kn 
At this moment the solution can be written 

N oo 

<f>(x,y) = -20 o ^ 



(19) 
(20) 



n=l /=— oo 
N 



h (v - Vn) - HD]' 



~ 2 ^°j^ E cosech2 (v ± X + C n)\ 



n=l 



14 



Regarding the determination of the constants c p , p = 1, N we note that 
the condition (J2Uj) only says that the differences between the values must be a 
multiple of m. However, there is the additional constraint that the constants 
cannot vanish. This is because the function would present singularities in 
the real space variable. On the other hand, the singularities must be placed 
symmetrically along the imaginary axis. These makes three restriction to 
any choice of the constants of integration 

1. the initial constraint, Eq.(|17|) which has been trasformed into Eq.(|20|) 

2. the restriction c p ^ 

3. the symmetrical positions along the imaginary axis, in order to have 
real solutions. 

This leads to the following choice 

in 

7Cfc = km + — 

where k = —N, ...,N. 

We will assume in the following that the number of poles is infinite, 

N -> oo 

However this will be discussed in more detail below. 

The problem of ennumeration of poles with ±x is so suppressed but we 
will have to make all steps separately for the two ± families. The fact that 
we can choose only two families is related to the impossibility to satisfy the 
constraints in the case when we would choose arbitrary combinations of the 
sign of x's in (y p - y n ). 

The constraints become two infinite systems of equations implying only 
Cfc, the constants. Then the solution can be written 

<t>{x,y) = -20 o — J^cosech 2 

71=1 
9 OO 

-20 O ^2"$^ COSeCh2 

71=1 

We can add to this expression the constraint equation, multiplied with (20 o ) 



7T7 



(y + x + c ri 



+ 



— (y - x + c' n ) 
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and written as an identity with zero. 

(f)(x,y) = 

9 OO 

~ 2 ^Jy2 1^ cosech 77 (v + x + c n) 



n=l 



+ 



+ ( 2< ^0) ^ E C ° SeCh2 ^ " C «)] " (^O) ( ^2 " 3 J 

n=l ^ ' 

o oo 

~ 2 ^~m 2^ cosech T> (y - x + c «) 



n=l 



2 00 / 2 i \ 

+ (20o) £j E cosech2 (s - <Q] - m (^-- 3 ) 

n=l ^ ' 

If the constants are chosen as suggested before then we have to exclude from 
the sum the term corresponding to n — 0. 



<f>(x,y) = 



-200^ E { cosech 



7T7 

— ( y + a; + — 



in nm 



-(20o 



71— — OO 

^ ^ o, ^ 12 

3^ " 3 j " 2< ^ C0Sech 



27 7 



7T7 



— cosech 



7T 

—nm 



D 



ITX 



— \y + x + — 



27 



-200^2" E | COSech2 



n=— oo 
n^O 

2 



27T n7T2 



7T7 

— I 7/ - X + — + 



(20 o ) 



7T 



3D 2 



2 0oT^ cosecl1 



27 7 



7T7 



cosech^ 



7T 

—nm 



17T 



For a reason that will become clear later we can only chose the positive sign 
of in/2. 
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This can further be written 

(f)(x,y) = 



-20o 



7T 



( 



\ 



+ ^ < cosech 



n=— oo 



7T7 

~D 



m nm 
y + x + — + 

2 7 7 



— cosech^ 



7T 

D 



nni 



7T 



-20 o -—cosech 2 



7T7 



y + X + 



27T 

27 



-20o^ 



7T 



- + < cosech 



V 



n=— oo 



7T 



-20 o -— cosech 2 
+3* 



7T 7 



D 



D 



in \ 



in nm 



"FT \ y - x + 7T + 



2 7 7 



cosech^ 



7T 

—nm 



\ y - x + 7T ) 



2 1 J 



Now we have to recall the identity for the doubly periodic elliptic Weier- 
strass function a 



1 . 2 f nw\ 

3 + C0S6Ch J + 



(21) 



+ < cosech 



n=— oo 



7T 



— (w + nLi) 

^2 



cosech 



= PW 



L2 

7T 



Where L\ and L 2 are respectively the period on the real axis and the period 
on the imaginary axis of the argument of p: 2u)\ = L\ and 2u>2 = iL 2 . We 
can identify 

L\ = ni 
L 2 = D 

and the two variables 



%7T 

w = i(y + x) + — 



(22) 



or w = 7 (y — x) + 
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At this moment an important comment should be done. As we have seen 
the solution of the constraint equation can be obtained on the basis of the 
comparison with the identity Eq. (jl8j) . We observe that the upper limit of 
the summation in ()18)) is interpreted as the number of poles, 

8-1 = N 

On the other hand we were led to identify 

D 

IT 

which means that the number of poles is the integer part 



Since we consider an infinite number of poles, N — > oo, this means that D is 
infinite. This is necessary since we later use the identity (121)1 where the sum- 
mation is extended over an infinite number of integer values n. Since however 
we find that L 2 = D, this implies that the period L 2 of the Weierstrass func- 
tion p is infinite. Certainly this is not acceptable in this particular approach, 
(but in general it is meaningful) and this shows that the identification using 
Eq. (j21j) can only be an approximation. This also means that the solution 
can only be approximative. However, if we decide to use this approximative 
identification, we have to understand in what consists this approximation 
and where we can expect to intervene the error we have introduced by that. 

We first note that the number of terms which is required in (J2"T|) to obtain 
a good approximation of the Weierstrass function is not necessarly large, for 
a significant area in the perodicity parallelogram of the complex argument. 
This means that actually D can be taken finite and in this case the use of 
Eq.(|18|) becomes legitimate. What is the number of poles, N and accord- 
ingly the length of periodicity Li = D in a reasonable approximation? From 
numerical experience it may be accepted that about five terms in the sum- 
mation determining the Weierstrass function gives a reasonable result. Then 
the number of poles retained in the sums (i.e. s — 1 = N) can be a few units 
and this also means that D ~few units. 

In the following we will write the solution as being expressed in terms of 
the Weierstrass function, but we have to remember that it is the result of an 
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approximation. It can be written 



(p(x,y) 



or 



" 2 0o^P 

TT 2 

- 2 0o^P 



Z7T 

7 (x + y) + — 



j(x-y) + 



ITY 



0, y) = -0o - 2 0o \ P 



l{x + y) + 



ITi 



+ P 



l{x-y) + 



(23) 



This expression represents in plane a system of rectangular cells. 

We will prove later that a single family of poles corresponding to one of 
the two possibilities ±x in Eq.([16p , which generates either the first or the 
second term in the Eq. (j23|) can provide an exact solution to the equation. 
This justifies the calculations presented in this section, since after them we 
are led directly to the form of the solution. 



7 The role of the elliptic function 



7.1 Basic information on the Weierstrass elliptic func- 
tion (lemniscate type) 

The most important consequence of the calculations that have been presented 
is the generation of a periodic solution and the confirmation, by the method 
of motion of poles, of the class of quasi-one-dimensional solutions. In the 
absence of a systematic integration procedure (like Inverse Scattering Trans- 
form) we can be at least sure that this class represents a significant part of 
the space of solutions. 

Even if Eq. fl23J) is the result of an approximation, it clearly shows that the 
solution should be searched in a form of a double periodic elliptic Weierstrass 
function. This function has appeared in our calculations independently of any 
attempt to integrate a one- dimensional version of the equation. 

We note that 



8P_ 

dx' 



-P 



l(x + y) + 



67 ( p 
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l(x + y) + 



ITi 
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Here g 2 is the second coefficient in the standard expression for the definition 
of the Weierstrass function 



p 1 (y) = u 



dt 



v/4t 3 - g 2 t - g 3 



In the following we will use this information to derive an exact solution 
of the Petviashvilli equation 



a<p - (3(j) 2 (24) 
We start by a linear substitution 

<j> (x, y) = sip (x, y) + — (25) 
(where s is a constant) which transform the equation into 

A "'=-* 2+ £H A (S) (26) 

We have allowed formally the space variation of the coefficients. For an easier 
reference to the properties of the Weierstrass function, we make a change of 
variables 

x — > x c = ix (27) 

y -> y c = iy 

which gives 

a 2 1 / a \ 
(x c , y c ) = (x c , yc)~^ + - A (-) (28) 

We look for a solution having the dependence on the coordinates x and 
y mediated by a new function, which we denote by u (x c , y c ) 

ip = ip(u) (29) 

with 

u = u (x c , y c ) (30) 
and now we proceed to express the equation (f2*Kj) using this form. 

dip dip du 
dx c du dx c 

d 2 ip d 2 ip ( du \ 2 dip d 2 u 



dx 2 du 2 \ dx c J du dx 2 
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and analogous for y. 

A c 4> (x c , y c 



d 2 i) 
du 2 

dip 



OU \ ( Oil 



+ 



du 



dx c 
[A c u] 



+ \dy c 



Suppose we find a function u (x c , y c ) verifying the conditions 





du 

dx n 



A c u 

2 



du\ 



where q is a constant. In this case the equation could be written 



d 2 ip sB , 9 a 2 1 . / a 

— - = — -ip 2 1 A — 

du 2 q Asj3q qs \23 



and this form has a known solution, from the identifications 



a 



s(3_ 
Q 

- —A (± 



6 

92 

2 



(31) 



(32) 



(33) 



As(3q qs \28 / 

This is because we have the known relationship for the Weierstrass function 



d 2 p(u) 2 g 2 
^u 2 ~ = ^ {U) 'J 

which has exactly the same form as (}3*2*|) . We find from Eq. ()33|) 

s/3 



(34) 



(35) 



92 



3a< 



6 . / a 
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(36) 
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8 Exact periodic solution of the Flierl-Petviashvilli 
equation 

We now turn to the solutions of the constraint equations defining u. It is 
clear that we may chose 

u (x c , y c ) = ay c + bx c + r (37) 
where a, b and r can be complex. The following condition results 

2 , 1.2 S @ 

a +b = — 
o 

In addition, a choice of a and b real numbers, which makes the first two 
terms in (j3*7|) purely imaginary, should be coroborated with the suggestion 
from Eq. (j22|) where the poles resulted shifted with a symmetric quantity 
with respect to the real axis. This time, due to the change of variables (j2*7|) 
everything is rotated. This yields the choice for r as half of the period (2a;) 
on the real axis 

(2^) 
2 

and the final form of the solution to the equation (j2H) 

(*, y) = - + s P ^ay + ib x + u% 2 = — 2 + —A ^- 



The second Weierstrass coefficient gs is left unspecified but it must be con- 
stant. 

We conclude that an exact solution to the Petviashvilli equation with 
constant coefficients a and f3 is 



a ( i \ 3a 2 \ 

< ' y) = ^ + s p y m v + lhx + = 7-^2 J ( :]> 



with the condition 

a 2 + b 2 = S 4 (39) 
o 

Now we understand the nature of the choice which is implicitely done 
when we discuss the solutions consisting of arbitrarly rotated ridges issued 
by translating one- dimensional profiles. It means to remain in the class of 
functions u whose Laplacean is zero and the gradient is constant, according 
to the conditions (|31|). Or, it can easily be seen that the only functions that 
verify these conditions have the form of linear combinations of the variables 
x and y. 
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9 The physical relevance of this periodic so- 
lution 



As a mathematical result, Eq.(|38J) is the exact response to the problem of 
solving the Eq.(j3J). This is valid for any sign of a and (3 since the function 
p is doubly periodic. The physical problem, as usual, is more complicated. 
The coefficients a and (3 which we have assumed constants actually have a 
certain space variation. We think however that it is worth examining the 
consequences of this solution for the nonlinear plasma models. 

In order to discuss possible physical applications of this solution we have 
to handle easily its numerical form. Although we have assumed that the two 
parameters a and (3 are constants, we will make estimations on the base of 
their physical origin. For this we remind that the basic physical constants 
are (Lakhin et al. J5], Spatscheck [Oj, Horton and Hasegawa pQ) 

a = 4fl--l (40) 



P 2 s V U 

T e 8 ( 1 \ e d ( 1 



2v?eBlp 2 s dx \L n J 2??vu 2 dx \L 

The coefficient a has the dimension (length)" 2 as it should. We introduce 
the following normalization for the potential 



T 

which only affects the coefficient (3 



d ( 1 \ T e 



2miU 2 dx \L n 
T P d f 1 



2uiiU 2 dx \ L 

d d ( l 



2u 2 dx \L 



At this moment the units are 



in the first term: A is measured in (m 2 ); the potential is adimension- 
alised, it is — > e<p/T e for which we have the order of magnitude 



n 
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• the second term: the coefficient a is in (to) 2 ; the potential is adimen- 
sionalised. 

• the third term: the coefficient (3 is in (m)~ 2 ; the potential is adimen- 
sionalised. 

With these values coefficients we have to calculate 

a 2 + 6 2 = f (41) 



92 = 7^ + ( - a ) (42) 



3a 2 6 ( a 

The latter parameter is connected with the periodicity lengths L\ and 1L2 or 
half-periods (u;,u/) of the Weierstrass function |T£j . 



. dt f°° dt 

u = / = = / . 43 

S1 \M 3 - <? 2 t - 03 4 y/4 (t - ei) (t - e 2 ) (t - e 3 ) 



C2 



e3 



y/4 (t - ei) (t - e 2 ) (t - e 3 ) 



A' 



The other period is 



= / ^ (44) 

A' 



In these formulas 



e = (45) 

ei - e 3 

is the modulus of the Jacobian elliptic functions and integrals. The similar 
quantity 

k' = VI - k 2 (46) 
is the complimentary modulus. The quantity 

K(k) = K = F(£,k) (47) 

de 

y/l - k 2 sin 2 9 
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is the complete elliptic integral of the first kind. And analogously 



K' = K (A/) 



(48) 



The order assumed is 



e 3 < e 2 < ei 



As will become later more clear, our problem has led to a definition of the 
Weierstrass function with only one coefficient , g 2 , fixed by the conditions, 
while gz cannot be made precise. This is because the equation we have used 
is the differential equation for the second derivative, p" (u) and this equation 
only depends on g 2 . 

The procedure for obtaining numerical results is : 

• assume physical parameters, density, temperature, etc. and calculate 

p s , v«, u, K n = L~\ dn n /dx; 

• calculate the coefficients of the original equation a and (3 from Eqs. (j40|) : 

• calculate g 2 from Eq.(jHJ) and chose the value of g 3 . Find ei, e 2 , e 3 ; 

• calculate the half-periods u and a/ on the real and imaginary axis from 
Eqs.flHl) and flSJ); 

• chose values for a and b such as to verify the Eq.(|41jl: 

• define the space region (y,x) i.e. (poloidal, radial) with y and x mea- 
sured in some typical Larmor radius p s0 ; 

• compute the complex variable: iay + ibx and scale to 2a/; take the real 
part of the argument as half the period on the real axis, in order to 
have real solution. 

• calculate the Weierstrass function of the complex argument, and asso- 
ciate it with the point (x, y), then calculate the solution by multiplying 
with s and adding the constant a/ (2/3); 



9.1 



An example 



We will assume aproximate values 



p s ~ 1CT 3 (m) 




0.5 
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Then physical values for these quantities are 



a 5 x 0.5 = 5 x 10 5 (m~ 2 ) 

no- 3 ) 2 v ; 



We take 



2 

~ 10 6 



and 

d_ 

dx 

then 



,8 ~ 10 7 (m" 2 ) 

We see that we can multiply all terms with p s ^ ~ 10~ 6 (m~ 2 ). In this 
way all distances will be expressed in units of Larmor radius p s0 . It results 

a — > a = 0.5 
/3 -> /3 = 10 

The space variation of the physical coefficients a and f3 are only in the 
radial x direction and is in general weak. The contribution of this part to 
the estimated value of g% is approximately one tenth from the first part (see 
also Appendix A) 

We chose a factor of scale s for the final amplitude of the part coming 
from p in the potential perturbation. This is a parameter that will result in 
general from the physical initial conditions, together with g 3 . We take 

s = 0.05 

For the following calculations, we leave g$ a free parameter. We use a 
computer code that calculates the value of the Weierstrass function for any 
complex argument, by reducing everything at a fundamental paralleogram 
with sides equal with 1 on the real as well as on the imaginary axis. This 
means that the real part of the argument must be scaled with (2a;) and the 
imaginary part of the argument must be scaled with (2a;'). After calculating 
the argument of the Weierstrass function, the value to be inserted in the 
subroutine is 

iay + ibx (1u)\ 1 



2a/ V 2 / 2uJ 
With these values the equations defining the parameters in our solution be- 



comes 

,2 , .2 



a z + b z = 0.09 
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We take 
and find 



3a 2 6 / a 



^3 = 



e 1 = 0.866 (49) 

e 2 = 

e 3 = -0.866 

k = 0.707 (50) 
k' = 0.707 

We find that the half-periods of the Weierstrass function are 

uj = 1.40879 (51) 

J = 1.40879 

We chose 

a = 0.0912 (52) 

b = 0.2738 

and we have to scale on the real axis with 

2co = 2.817 (53) 

and on the imaginary axis 

2u' = 2.817 (54) 
This gives a first estimation of the width of the layer along the x direction: 

2uj' 

Sx ^ _ ^ 10 (55) 

which means about 1 cm. We calculate the profile of the streamfunction 
4>(y,x) on a poloidal-radial domain of extension (40p s x 40p s ). 

For example, we show the structure of the solution in the two plots. 

It is interesting to examine the structure of the flow field induced by the 
potential as a function of the position along the minor radius in tokamak. We 
first perform a numerical simulation using a one dimensional transport code 
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Figure 1: The streamfunction <\> solution of the Petviashvilli equation for the 
ratio \ = 0.9. 
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-10 10 

y (poloidal) in units of p 



Figure 2: The same as Figure 1, in contour plot. 
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|25j . Details are given in Appendix B. The physical data from simulation 
consists of 

T e (r) , Ti (r) , n (r) , p s (r) , c s (r) , L n (r) , L Te (r) , (r) 

which are transfered to a code that calculates the parameters a (r) , /3 (r), 
a (r), 6 (r), ^2 (?") and, takes a fixed value for #3. Then the half-periods uo (r) 
and u/ (r) are calculated. These parameters are represented in the figures 
below. We have studied the effect of the variation of the second Weierstrass 
parameter, g 3 . When it is taken equal to zero, we find that the width of the 
layer of poloidal flow, Eq. (J33J) is proportional to the local Larmor radius, 

6x = 3.924K (k') (1 - v*/ U y 1/2 p s 

(with the constants of Eqs. (|48j) and (|5()jl ). The variations of the width of 
a layer Sx/p s and of the maximum perturbation e<p/T e with g 3 and v*/u are 
shown in FigsEl and [7| On this graphs only the part where the discriminant 
of the Weierstrass function is positive A = g\ — TJg\ > are shown, since 
there the roots 61,2,3 are all real (these points are also indicated in projection 
on the plane). The characteristics of an initial physical field that can evolve 
toward this layered geometry of flow can be inferred from these graphs by 
eliminating the common variable g 3 . 

10 Conclusion 

It has been known from the experiments on Rayleigh-Benard convection that 
at higher Reynolds number there is a second bifurcation (the first being from 
purely conducting to convective rolls phase) where the pattern of the fluid 
flow exhibits tilted structures. Close to the upper and to the lower boundaries 
there is a deformation of the rolls, the "wind" as has been designated, which 
significantly enhances the Reynolds stress and favorizes second instability 
and eventually global displacement of the fluid along the boundaries. This 
finding, supported independently by suggestive results from the numerical 
simulations of the ITG instability led to the idea that a tilting instability 
may favorise, by an increase of the Reynolds stress, the generation of zonal 
flow in plasma. The solutions we have found is a new factor that should be 
included in this physical description. 

Certainly, the tilt of the radially elongated eddies specific to the ITG 
potential pattern can be favorable to the spontaneous generation of sheared 
poloidal flow. On the other hand, the solution we have found has a strong 
similarity to the tilted cells and zonal flows. It may be possible that after a 
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Figure 3: Physical parameters and the derived values a , (3, a and b as 
function of r. 
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Figure 5: Width of the layer of flow and amplitude induced by the pertur- 
bation. 
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Figure 6: Variation of the amplitude e0/T e with g 3 and v*/u. 



60- 




-2 0.4 

Figure 7: Variation of the width 8x of a layer with g 3 and t>*/w. 
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certain degree of tilting is attained, the plasma evolves spontaneously to a 
solution of the type described above. This solution intrinsically consists of 
layers of sheared flow. What is more important is that this solution has all 
the usual attributes of an exactly integrable structure: it is more robust and 
may represent an attractor. It can be considered that relating the generation 
of the zonal flow to the process of evolution of a system to a robust attractor 
is a useful approach to be further examined. 

In this work we have restricted to the nonlinear dynamics of the 2D 
structures. As a further extension of the nonlinear solution of electrostatic 
drift waves described by the Hasegawa-Mima equation, the electromagnetic 
3D vortical motion has also been studied 26J, [27J. In particular it has been 
shown that a 2D vortical motion in the (r, #)-plane propagates along the 
magnetic line of force. Nonlinear coupling of drift-shear Alfven waves may 
lead to a novel nonlinear solution. Recently this topic has been extended and 
the formulation has been updated in Ref. [28J. These are however out of the 
scope of the present work and will be a subject of future research. 
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Appendix A : Discussion on the weak space 
variation of the parameters 

In order to obtain the exact solution we had to assume that the coefficients 
are constant. Later, we have considered a weak variation of the coefficients 
with the coordinate x, induced by the presence of the physical parameters. 
We can assume that this is a reasonable representation of the real situation 
only if we have an adiabatic variation of gi with the physical parameters. 

The numerical study of the dependence of the elliptic function on the 
physical parameters shows that g% is indeed slowly varying and that the 
correction due to the second term in the expression of g% can be neglected. 

One can see that the weaker restriction can be formulated as follows: g% 
is independent on« = iay + ibx + ui, which means that #2 is constant along 
the lines ay + bx = const and can only have a slow parametric variation 
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in the direction transversal to this family of lines. Taking into account the 
expression of g 2 and that a and j3 has variation mainly along the radial 
direction, we can reformulate the restriction by requiring that g 2 has no 
variation along lines quasi-parallel to the poloidal dierction, and can only 
have a weak dependence perpendicular on these lines. This explains our 
choice of combination of a and b in the numerical study presented above. 

We will make a scaling transformation to reduce the Weierstrass function 
to constant coefficient g 2 . For this we recall the formula for homogeneity of 
the Weierstrass function 



p(u;g 2 ,g 3 ) 



We take 



H 2 P 



33 = 



92 93 

fi 4 fl 6 



(A.l) 



1/4 
V = 92 



3a 2 



s 2 (3 2 s 2 (3 \/3 



a 



1/4 



and the formula becomes 

p («; 92, gs) = g 2 2 p (g 2 A u- 1, o) 



or 



(A.2) 



a 



gl /2 sp (g 2 4 (iay + ibx + u) ; 1, o) (A.3) 



This formula may serve for a numerical investigation of the space varia- 
tion and the estimation of the error in using the solution based on constant 
coefficients. 

We ask the stronger condition, that g 2 is simply a constant 



92 



3a 2 



s 2 f3 2 s 2 (3 \P 



const 



This is essentially a differential equation which strongly constrain the space 
dependence of the physical parameters. Taking for example const= p, and 
remembering that the physical variation is on x, we have 



d 2 ( a 
d^ 2 \J3 



a 
2? 



(3s 2 p 
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Now we assume that the radial third derivative of the density is small and 
take then j3 constant. The equation for a is 



a 



Qa z - 2(3 2 s 2 p 



(A.4) 



after a dividing the variable x with vl2. This is again the Weierstrass 
equation and the solution is p (x) where x is measured from the surface 
where the periods are calculated. Looking for non-periodic solutions we find 



a (x) = 1 



u 



(3 2 s 2 p 



1/2 



3 (/?Vp/3) 



1/2 



cosh^ 



\/3(/?Vp/3) 



1/4 



(A.5) 



We see that the behaviour of the right hand side in Eq. (|A.5|) is approximately 

linear in the variable (/3 2 s 2 p/3) 1 ^ 2 . Then we can expect a condition of the 
form 

v* „ (92\ x l 2 



u V 3 



or 



L n 2udx\L n J \3/ 



Here all distances, x and L n have been normalized at a typical Larmor radius, 
p s0 . This can be reduced to a condition on the density variation with the 
minor radius 



x<>"- 

dr 



3 / 2u 



0.05 



The conclusion is that the density gradient length has a variation of the type 



exp -0.05- 



x 



PsO 

on an interval < x < 100p s0 in the region where this solution exists. The 
fast variation of the density is favorable to the validity of this solution. 



Appendix B : Numerical simulation for toka- 
mak plasma parameter's profiles 

The code solves the balance equations for energy, density and fields. The 
variables are: the electron and ion temperatures T e , Tf, the current density 
j, the poloidal magnetic field B e , the toroidal electric field E v , the electron 
density n e , the radial pinch particle velocity V r . The following equations 
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are discretized on a one-dimensional space mesh (on the small radius) and 
evloved in time by a semi-implicit scheme. 



3 d_ 
2dt 





+ n e VT e ))+Ej 



add 



3d_ 
2dt 



(riiTi) 




1 1 d 



(rB e ) 



fi r dr 



8B e _ BE 
dt dr 

E = m 

dn e _ ld_ ( dn e 
dt r dr \ dr 



Neutral atoms as well as impurities (Carbon, Oxygen, Iron, Wolfram, Molib- 
den) are considered. The transport coefficients are those of the Merejkhin- 
Mukhovatov model and the specific parameters of the tokamak correspond 
to the Tore-Supra device. The run is extended over 20 seconds and the sta- 
tionarity is reached within few seconds. Then we extract radial dependent 
plasma parameter from a time slice at about 13.5 seconds. 
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